unit frmMainUnit;
{$POINTERMATH ON}

interface

uses
  Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants,
  System.Classes, Vcl.Graphics, Vcl.Controls, Vcl.Forms, Vcl.Dialogs,
  Vcl.StdCtrls, Math, VCLTee.TeEngine, Vcl.ExtCtrls, VCLTee.TeeProcs,
  VCLTee.Chart, VCLTee.Series;

type
  TfrmMain = class(TForm)
    Button3: TButton;
    Memo1: TMemo;
    Button4: TButton;
    procedure Button3Click(Sender: TObject);
    procedure Button4Click(Sender: TObject);
  private
    procedure CalcTables(AResolution, ABitsA: Byte; AGenSource: Boolean;
                         out ASize: Integer; out AMaxError: Double;
                         out S: string);
  end;

var
  frmMain: TfrmMain;

implementation

{$R *.dfm}

procedure TfrmMain.CalcTables(AResolution, ABitsA: Byte; AGenSource: Boolean;
  out ASize: Integer; out AMaxError: Double; out S: string);
var
  A_BITS: UInt64;
  B_BITS: UInt64;
  A_SIZE: UInt64;
  B_SIZE: UInt64;
  DBL_SZ: UInt64;
type
  TItem = record
    A: Double;
    B: Double;
  end;
  PTabB = ^TTabB;
  TTabB = record
    Next:       PTabB;
    StartIndex: Integer;
    Hash:       Integer;
    Data:       array of TItem;
  end;

  PInts = ^TInts;
  TInts = record
    I1, I2: Integer;
  end;

  function _Sigmoid(X: Double): Double;
  begin
    if x <= -2 then
      Result := -1
    else if X >= 2 then
      Result := 1
    else
      Result := 2 / (1 + exp(- 4.9 * X)) - 1;
  end;

  function _IsEqual(T1, T2: PTabB): Boolean;
  var
    j: Integer;
  begin
    for j := 0 to B_SIZE - 1 do
      if (T1.Data[j].A <> T2.Data[j].A) or
         (T1.Data[j].B <> T2.Data[j].B) then
        Exit(False);
    Result := True;
  end;

  function _Find(T, ATabs: PTabB): PTabB;
  var
    cur: PTabB;
  begin
    cur := ATabs;

    while cur <> nil do
    begin
      if (cur.Hash = T.Hash) and _IsEqual(cur, T) then
        Exit(cur);
      cur := cur.Next;
    end;

    Result := nil;
  end;

var
  tabA:    array of PTabB;
  tabsB:   PTabB;
  tabBcnt: Integer;
  i, j:    Integer;
  idx:     Integer;
  cur:     PTabB;
  x1, x2:  Double;
  y1, y2:  Double;
  fnd:     PTabB;
  xe:      Double;
  err:     Double;
  eitm:    TItem;
begin
  A_BITS := ABitsA;
  B_BITS := ARESOLUTION - A_BITS;
  A_SIZE := 1 shl A_BITS;
  B_SIZE := 1 shl B_BITS;
  DBL_SZ := SizeOf(Double) * 8;

  { Generate tables }

  SetLength(tabA, A_SIZE);
  tabsB     := nil;
  tabBcnt   := 0;

  for i := 0 to A_SIZE - 1 do
  begin
    New(cur);
    SetLength(cur.Data, B_SIZE);
    cur.Hash := 0;

    for j := 0 to B_SIZE -1 do
    begin
      PUInt64(@x1)^ := (UInt64(i) shl (DBL_SZ - A_BITS)) or
                       (UInt64(j) shl (DBL_SZ - ARESOLUTION));
      PUInt64(@x2)^ := PUInt64(@x1)^ or ((UInt64(1) shl (DBL_SZ - ARESOLUTION)) - 1);

      // Two points line equation:
      //   Y = (y2 - y1) / (x2 - x1) * X + (x2 * y1 - x1 * y2) / (x2 - x1)

      if not IsNan(x1) and not IsNan(x2) then
      begin
        y1            := _Sigmoid(x1);
        y2            := _Sigmoid(x2);
        cur.Data[j].A := (y2 - y1) / (x2 - x1);
        cur.Data[j].B := (x2 * y1 - x1 * y2) / (x2 - x1);
      end
      else
      begin
        cur.Data[j].A := 0;
        cur.Data[j].B := 0;
      end;

      cur.Hash := cur.Hash xor PInts(@cur.Data[j].A).I1 xor
                  PInts(@cur.Data[j].A).I2 xor PInts(@cur.Data[j].B).I1 xor
                  PInts(@cur.Data[j].B).I2;
    end;

    fnd := _Find(cur, tabsB);
    if fnd <> nil then
    begin
      Dispose(cur);
      cur := fnd;
    end
    else
    begin
      cur.Next := tabsB;
      tabsB    := cur;
      Inc(tabBcnt);
    end;

    tabA[i] := cur;
  end;

  { Calc error and size }

  AMaxError := 0;
  xe        := -20;

  while xe <= 20 do
  begin
    eitm := tabA[(PUInt64(@xe)^ shr (DBL_SZ - A_BITS))].Data[
                    (PUInt64(@xe)^  shr (DBL_SZ - ARESOLUTION)) and (B_SIZE - 1)];
    err  := Abs(_Sigmoid(xe) - (eitm.A * xe + eitm.B));

    if err > AMaxError then
      AMaxError := err;
    xe := xe + 0.01;
  end;

  ASize := A_SIZE * SizeOf(Pointer) + tabBcnt * B_SIZE * 2 * SizeOf(Double);

  { Generate source }

  if AGenSource then
  begin
    S := S + '// Piecewise Sigmoid function approximation: F(X) = 1 / (1 + Exp(-X)).'#13#10;
    S := S + '// ==================================================================='#13#10;
    S := S + '//   Precision:          ' + IntToStr(AResolution) + ' bits.'#13#10;
    S := S + '//   Max error detected: ' + FloatToStr(AMaxError) + '.'#13#10;
    S := S + '//   Overall data size:  ' + IntToStr(ASize) + ' bytes.'#13#10;

    S := S + #13#10;
    S := S + 'SHIFT_A = ' + IntToStr(DBL_SZ - A_BITS - SizeOf(Integer) * 8) + ';'#13#10;
    S := S + 'SHIFT_B = ' + IntToStr(DBL_SZ - ARESOLUTION - SizeOf(Integer) * 8) + ';'#13#10;
    S := S + 'MASK_B  = $' + IntToHex(B_SIZE - 1, 8) + ';'#13#10;


    S := S + #13#10;
    S := S + 'SUBTBL: array[0..' + IntToStr(tabBcnt * B_SIZE - 1) +
         '] of TLine = ('#13#10;

    idx := 0;
    cur := tabsB;

    while cur <> nil do
    begin
      cur.StartIndex := idx;
      FormatSettings.DecimalSeparator := '.';

      for i := 0 to B_SIZE - 1 do
      begin
        S := S + '  (A: ' +  FloatToStrF(cur.Data[i].A, ffExponent, 100, 3) +
             '; B: ' + FloatToStrF(cur.Data[i].B, ffExponent, 100, 3) + ')';
        if (i <> B_SIZE - 1) or (cur.Next <> nil) then
          S := S + ',';
        S := S + #13#10;

        Inc(idx);
      end;

      cur := cur.Next;
    end;

    S := S + ');'#13#10;

    S := S + #13#10;
    S := S + 'DATA: array[0..' + IntToStr(A_SIZE - 1) +
         '] of PLine = ('#13#10;

    for i := 0 to A_SIZE - 1 do
    begin
      S := S + '  @SUBTBL[' + IntToStr(tabA[i].StartIndex) + ']';
      if i <> A_SIZE - 1 then
          S := S + ',';
        S := S + #13#10;
    end;

    S := S + ');'#13#10;
  end;
end;

procedure TfrmMain.Button3Click(Sender: TObject);
var
  bestsz, sz: Integer;
  maxerr:     Double;
  p, btsa:    Integer;
  bestbtsa:   Integer;
  s: string;
begin
  Memo1.Lines.Clear;

  for p := 3 to 16 do
  begin
    bestsz   := MaxInt;
    bestbtsa := -1;

    for btsa := 1 to p - 1 do
    begin
      CalcTables(p, btsa, False, sz, maxerr, s);
      if sz < bestsz then
      begin
        bestsz   := sz;
        bestbtsa := btsa;
      end;
    end;

    Memo1.Lines.Add('=== Precision: ' + IntToStr(p) + ' ===');
    Memo1.Lines.Add('  BitsA: ' + IntToStr(bestbtsa));
    Memo1.Lines.Add('  Tables Size: ' + IntToStr(bestsz));
    Memo1.Lines.Add('  Max Error: ' + FloatToStr(maxerr));
    Memo1.Update;
  end;
end;

procedure TfrmMain.Button4Click(Sender: TObject);
const
  PRECISION = 14;
  A_BITS    = 9;
var
  sz:     Integer;
  maxerr: Double;
  s:      string;
begin
  s := '';
  CalcTables(PRECISION, A_BITS, True, sz, maxerr, s);
  Memo1.Lines.Text := s;
end;

end.
